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Abstract 

We have performed a direct numerical simulation of dilute turbulent particulate flow in a 
vertical plane channel, fully resolving the phase interfaces. The flow conditions are the same 
as those in the main case of "Uhlmann, M., Phys. Fluids, vol. 20, 2008, 053305", with the 
exception of the computational domain length which has been doubled in the present study. 
The statistics of flow and particle motion are not significantly altered by the elongation of 
the domain. The large-scale columnar-like structures which had previously been identified 
do persist and they are still only marginally decorrelated in the prolonged domain. Voronoi 
analysis of the spatial particle distribution shows that the state of the dispersed phase can 
be characterized as slightly more ordered than random tending towards a homogeneous spa- 
tial distribution. It is also found that the p.d.f.'s of Lagrangian particle accelerations for 
wall-normal and spanwise directions follow a lognormal distribution as observed in previous 
experiments of homogeneous flows. The streamwise component deviates from this law pre- 
senting significant skewness. Finally, a statistical analysis of the flow in the near field around 
the particles reveals that particle wakes present two regions, a near wake where the velocity 
deficit decays as a;~^ and a far wake with a decay of approximately a;~^. 

Keywords: particulate flow, direct numerical simulation, interface resolution, turbulent chan- 
nel flow 

1 Introduction 

Fluid flow with suspended solid particles is encountered in a multitude of natural and industrial 
systems. Examples include the motion of sediment particles in rivers, fluidized beds and blood flow. 
Despite the great technological importance of these systems our understanding of the dynamics 
of fluid-particle interaction is still incomplete at the present date. Recently, however, significant 
progress has been made based on data provided by new experimental methods as well as numerical 
simulations. While most past investigations of numerical type have been performed in the context 
of the point-particle approach, it has now become possible to simulate the motion of a considerable 
number of finite-size particles including an accurate description of the surrounding flow field on 
the particle scale (Pan and Banerjee, 1997; Kajishima and Takiguchi, 2002; Ten Gate et al., 2004; 
Uhlmann, 2008; Lucci, Fcrrante and Elghobashi, 2010, 2011). Although the complexity of these 
particle-resolved simulations (in terms of Reynolds number, number of particles and computational 
domain size) is still limited, new insight into the physics of fluid-particle systems is beginning to 
emerge from such studies. 

Uhlmann (2008) has simulated turbulent flow in a vertically-oriented plane channel seeded 
with heavy spherical particles with a diameter corresponding to approximately 11 wall units at a 
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solid volume fraction of 0.4%. The pressure-driven upward flow (at constant flow rate) was found 
to be strongly modified due to the particle presence, with increased wall-shear stress and strongly 
enhanced turbulence intensity. The average relative flow, corresponding to a Reynolds number 
(based on particle diameter) of approximately 135, lead to the establishment of wakes behind 
individual particles. Additionally, the formation of very large-scale, streak-like flow structures 
(essentially spanning the entire box-size) , absent in corresponding single-phase flow, was observed. 
At the same time the dispersed phase did not exhibit any of the common signs of preferential 
concentration. 

In the present study we are revisiting the same flow conflguration of vertical particulate channel 
flow, expanding upon the previous analysis of Uhlmann (2008) by addressing several unanswered 
questions. First, we wish to determine the influence of the streamwise length of the computa- 
tional domain upon the largest flow scales. For this purpose we have performed new simulations 
analogous to the ones conducted by Uhlmann (2008), but with twice the value of the original 
streamwise period, while keeping all remaining parameters unchanged. 

Second, we intend to provide a more complete description of the turbulent fluid-particle interac- 
tion in vertical channel flow. To this end we have analyzed three aspects of the flow dynamics which 
had previously not been considered by Uhlmann (2008): Voronoi analysis of the spatial structure of 
the dispersed phase, analysis of particle acceleration statistics, and particle-conditioned averaging 
of the fluid flow field. 

Voronoi analysis is a relatively recent addition to the arsenal of tools for the description of 
particles suspended in fluids (Monchaux, Bourgoin and Cartellier, 2010). In the present flow 
configuration it turns out that this methodology provides a more sensitive measure of the particle 
phase geometry than previously employed criteria. 

The statistical properties of particle acceleration have received increasing attention in recent 
years (Toschi and Bodenschatz, 2009). Since particle acceleration is (up to particle mass) equiv- 
alent to the resulting forces acting upon the particles, its analysis can be instrumental in un- 
derstanding turbulence-particle interaction mechanisms. One application where the infiuence of 
turbulence upon particle acceleration statistics is believed to be of key importance is the growth 
of rain drops by collisions in atmospheric clouds (Warhaft, 2009). Modern experimental results on 
the acceleration of finite-size particles (Qureshi et al., 2007; Xu and Bodenschatz, 2008; Brown, 
Warhaft and Voth, 2009) have only started to emerge around the date of publication of the pre- 
cursor paper (Uhlmann, 2008). Therefore, such an analysis was not carried out therein. Here we 
present a statistical analysis of particle acceleration/hydrodynamic forces, relating the findings to 
available experimental results. 

Finally, the understanding of the interaction between solid particles and fiuid turbulence does 
not seem complete without a statistical analysis of the flow in the near-field around the particles. 
In order to investigate the characteristics of particle-induced wakes and with the aim to provide 
data which might be useful for the purpose of two-phase flow modelling, we have undertaken a 
study of particle-conditioned averaging of the flow field. Reference data for fixed particles swept 
by (essentially) homogeneous-isotropic flow (Bagchi and Balachandar, 2004; Amoura et al., 2010) 
as well as wall-bounded shear fiow (Wu and Faeth, 1994; Legendre, Merle and Magnaudet, 2006; 
Zeng, Balachandar and Najjar, 2010) is available and has been used for the purpose of comparison. 

2 Computational setup 
2.1 Numerical method 

The numerical method employed in the current simulations is identical to the one detailed in 
Uhlmann (2005) which was already used for the previous simulations of vertical particulate channel 
flow by Uhlmann (2008). The incompressible Navier-Stokes equations are solved by a fractional 
step approach with implicit treatment of the viscous terms (Crank-Nicolson) and a three-step 
Runge-Kutta scheme for the non-linear terms. The spatial discretization employs second-order 
central finite-differences on a staggered mesh. The no-slip condition at the surface of moving solid 
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Figure 1: Illustration of the computational domain, which is bi-periodic in the streamwise {x) and 
spanwise (z) directions, (a) shows the domain used in Uhlmann (2008), (b) the current domain 
which is a streamwise extension of the former (by a factor of two) . The red spheres indicate actual 
instantaneous particle positions. 

Reb Rer fj \g\h/ul St+ Sh h/D D+ 

2700 220.9 2.2077 12.108 15.5 0.83 20 11.25 0.0042 

Table 1: Physical parameters for particulate flow in a vertically oriented plane channel: bulk 
Reynolds number Reb (imposed quantity), friction- velocity based Reynolds number Rer (derived 
quantity), particle/fluid density ratio Pp/pf, gravitational parameter |g|/i/u^. Stokes numbers 
based upon bulk units Stb = TpUb/h (imposed) and wall units St'^ = TpU^/v (derived), length 
scale ratio h/D (imposed) and D+ (derived), and global solid volume fraction $s. The particle 
relaxation time scale was defined as Tp = ppD'^/{pf 18 ly). 



particles is imposed by means of a specially designed immersed boundary technique (Uhlmann, 
2005) . The motion of the particles is computed from the Newton equations for linear and angular 
motion of rigid bodies, driven by buoyancy, hydrodynamic forces/torque and contact forces (in case 
of collisions). Since the suspension under consideration is dilute, collisions are treated by a simple 
repulsive force mechanism (Glowinski et al., 1999) formulated such as to keep colliding particles 
from overlapping non-physically. The same treatment is applied to particle-wall encounters. It 
should be noted that the employed computational grid is uniform and isotropic. The chosen grid 
width Aa; = Ay = Az yields a particle resolution of Dj Ax — 12.8, a resolution of the channel 
half-width of hj Ax = 256 and Ax"*" = 0.67 in terms of wall units. Further information on our 
extensive validation tests and grid convergence can be found in Uhlmann (2005, 2008) and further 
references therein. 
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case Cl X Ny X Nz Np tobsUb/h 

Uhlmann (2008) 8h x 2h x Ah 2048 x 513 x 1024 4096 115 
present l&h x 2h x Ah 4096 x 513 x 1024 8192 90 

Table 2: Numerical parameters employed in the simulations: computational domain size fl, number 
of grid nodes Ni in the ith coordinate direction, number of particles Np, temporal observation 
interval tabs after discarding the initial transient. The grid spacing in all cases is fixed at Ax — 
h/256, corresponding to Nj^ ~ 515 Lagrangian force points distributed over the surface of each 
particle. 



(a) (6) 
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Figure 2: Two-point correlation functions of fluid phase velocity fluctuations for streamwise sepa- 
rations computed from a single snapshot of the flow field in a wall-parallel plane at a wall distance 
of y/h = 0.0938 (i.e. y+ ss 22). (a) streamwise velocity component Rmi, (b) spanwise velocity 
component R^w ■ The correlation immediately after periodic extension from the small-box simula- 
tion of Uhlmann (2008) is shown as a dashed line; the correlation at the beginning of the present 
sampling interval is indicated by a solid line. 



2.2 Flow configuration 

Figure 1 shows the flow configuration under consideration as well as the coordinate system. The 
plane channel is oriented vertically, x being the streamwise coordinate direction, y is the wall- 
normal (with the channel width equal to 2h) and z the spanwise direction. Fluid flow is directed 
upwards (in the positive x direction), driven by a streamwise pressure-gradient. The bulk velocity 
Ub is maintained at a constant value, such that the Reynolds number based upon the bulk velocity, 
Rcb — Ubh/v, is imposed (cf. table 1 for the values of the principal physical parameters). A large 
number {Np = 8192) of monodispersed, rigid, spherical particles is suspended in the flow. The 
nominal terminal velocity of the particles (computed from an equilibrium of buoyancy force and 
standard drag force, Clift, Grace and Weber, 1978) is set equal to the bulk velocity of the fluid 
phase. Consequently, the average particle settling velocity obtained in the actual simulation 
is roughly zero. The chosen density ratio Pp/pf = 2.2077 (where Pp, pf are the particle and 
fluid densities) is comparable to the case of glass particles in water. The particle diameter D, 
approximately equal to 11 wall units, is comparable to the cross-sectional scales of buffer layer 
flow structures. Finally, table 1 shows that the suspension is indeed dilute, with less than one half 
percent of solid volume fraction. 

As can be seen from table 2, the present simulation is performed in a computational domain 
which has twice the streamwise period as compared to Uhlmann (2008), while maintaining an 
identical small-scale resolution. The table also shows that an observation interval of approximately 
90 bulk time units (defined as Tb — h/ub) has been computed after discarding the initial transient. 
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Figure 3: (a) Mean velocity profiles of both phases. The fluid phase is shown as: , present; 

, Uhlmann (2008). The mean velocity of the particulate phase is shown as: •, present; o, 

Uhlmann (2008). (b) Apparent slip velocity between the phases: •, present; o, Uhlmann (2008). 
Fluid data shown in this figure is 'composite-averaged' (cf. equation 23) for the purpose of strict 
comparison; pure fluid averaging data for the present simulation is shown in Appendix B. 

which will be discussed next. 

The simulations were run on different supercomputing systems, typically using between 256 
and 1024 processor cores. The total number of CPU hours spent was of the order of 5 million. 

2.3 Simulation start-up and initial transient 

The current simulation was initialized with an exact periodic extension of a flow field taken at an 
instant towards the end of the simulation of Uhlmann (2008). For this purpose, fluid and particle 
data in the interval x G [0, 8h] was copied from the reference field. In order to obtain data in 
the interval x € [8h, 16h], the shift x = x + 8h was applied to the reference field. It should be 
emphasized that no explicit perturbations whatsoever were added to the initial field. 

Subsequently, the extended simulation was run while different quantities were monitored in 
order to determine whether the system has developed sufficiently such as to "forget" the initial 
state. A sensitive measure of independence from the initial condition is provided by two-point 
correlations of fluid data. Therefore, correlation functions of fluid velocity components as a func- 
tion of streamwise separations have been analyzed. As can be seen from figure 2, the field initially 
exhibits a periodicity over half the domain length, consistent with the fact that the simulation 
was started from an exact "copy" of a field taken from a run in a domain with half the stream- 
wise period. Over an interval of approximately 40 Tf, the artificial periodicity (with period 8h) 
gradually disappears and only the strict periodicity over the fundamental period of 16ft. remains. 
The transient interval was discarded and statistics have then been accumulated over a sampling 
interval of 90 Tb. The plots in figure 2 also show the correlation function computed from a snapshot 
at the beginning of the interval over which statistics have been accumulated, which clearly lacks 
any traces of artificial periodicity. 

2.4 Notation 

Before turning to the results, let us fix the basic notation followed throughout the present text. 
Velocity vectors and their components corresponding to the fluid and the particle phases are distin- 
guished by subscripts "f" and "p", respectively, as in u/ = {uf,Vf,Wf)'^ and Uj, = {up,Vp,Wp)'^ . 
Similarly, the vector of angular particle velocity is denoted as ujp = (wpi, ajp2, i^ps)^ and lin- 
ear particle acceleration as = (ftpi, ap2, ips)"^, while fluid acceleration is denoted as af = 
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Figure 4: Wall normal profiles of: (a) the mean solid volume fraction, (6) the mean value of the 
spanwise component of angular particle velocity. Symbols as in figure 3(5). 



{afi,af2,aft^)'^ . The equations of motion for an immersed sphere with density Pp can be written 
as: 

PpVpAp = Fh+Fb + Fc, (la) 
Ipibp = Th, (lb) 

where Vp = ttD^ /6 is the volume of a sphere with diameter D and Ip = ppD^n/GO is its moment 
of inertia, considering a homogeneous mass density. In equation (la) the resulting force has been 
separated into a hydrodynamic contribution Fh — r ■ nda — J^pnda {S being the surface of 
the sphere, n the outward pointing normal vector at the surface, t = pfV {dufi/dxj + dufj/dxi) 
the viscous stress tensor and p the hydrodynamic pressure), a relative buoyancy force — 
{pp — pf)gVp and a contribution Fc from solid-solid contact which is modelled as discussed in 
§2.1. The angular acceleration in (lb) only has a single contribution from viscous hydrodynamic 
stresses, viz. Th = x (r • n) da {r^ being the distance vector from the particle center), since 
our solid-solid collision treatment is limited to normal forces. 



3 Results 

3.1 Effect of streamwise elongated domain 

The aim of the present section is to examine the effect of the finite streamwise period imposed 
upon the fields in the simulation. For this purpose we will compare the present results with those 
obtained in a shorter domain as presented by Uhlniann (2008). 

Note that the fluid data shown in figures 3 and 5 are based upon 'composite' averaging, i.e. 
not distinguishing between velocity values of those grid nodes which are instantaneously located 
in the fluid phase and those inside the solid phase. Doing so allows for a strict comparison with 
results from the simulation of Uhlmann (2008), where this distinction was not made (cf. discussion 
in Appendix A of Uhlmann, 2008) . The specific fluid-phase statistics (including only actual grid 
nodes located inside the fluid domain, cf. definition in equation 22) corresponding to the fluid 
quantities shown in figures 3 and 5 are included in Appendix B for future reference. 

In the present flow it turns out that the average time interval between two collision events 
experienced by a particle is 3.9 Tf,, showing that indeed collisions are relatively infrequent. Fur- 
thermore, no particle-wall collisions have been observed during the simulation interval. 
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Figure 5: (a) R.m.s. of velocity fluctuations of both phases, with lines corresponding to the fluid 

phase ( , present; , Uhlmann (2008)), and symbols to the particulate phase (•, present; 

o, Uhlmann (2008)), as in figure 3(a). The color coding indicates streamwise (black), wall-normal 
(red), spanwise (blue) components, (b) Reynolds shear stress of fluid velocity fluctuations as well 
as corresponding velocity correlation of the particle motion. All quantities are normalized in bulk 
units. Fluid data shown in this figure is 'composite-averaged' (cf. equation 23) for the purpose of 
strict comparison; pure fluid averaging data for the present simulation is shown in Appendix B. 

3.1.1 Eulerian statistics of both phases 

Mean fluid and particle velocity profiles as well as particle concentration profiles and the mean 
spanwise component of the angular particle velocity are shown in figures 3 and 4. It can be 
seen that these average values are not too strongly affected by an increase of the box size from 
Lx/h = 8 to 16. The observed mostly small differences in these quantities can be attributed to the 
remaining statistical uncertainty (i.e. limited sampling of the largest ffow scales). In particular, 
the wall shear-stress averaged over the observation interval differs slightly, such that the time-and- 
wall-averaged friction Reynolds number Rcr takes a value of 220.9 (versus 224.4 in the simulation 
of Uhlmann, 2008) . Furthermore, the previously observed weak tendency towards a concave mean 
velocity profile (Uhlmann, 2008) is not confirmed by the present results (cf. figure 3a). It might 
therefore equally be attributed to limited sampling of the large scales in the previous simulation. 
The difference between the mean velocities of both phases (cf. figure 36) is termed the 'apparent 
slip velocity', which will be denoted by uiag — (up) — (uf) in the following. The corresponding 
Reynolds number {Reiag — \uiag\D/iy) measures Reiag ~ 132 in the bulk of the channel, while the 
value drops significantly for wall-distances smaller than y/h ^ 0.1. 

Concerning the mean value of the spanwise component of angular particle velocity shown in 
figure 4(6), it should be mentioned that the corresponding graph in Uhlmann (2008) erroneously 
shows the quantity {ujp^z)ub/h (i.e. the axis label of figure 14 in Uhlmann, 2008, is incorrect). 
The correct non-dimensional quantity {ujp z)h/uf, from that reference is shown in the present fig- 
ure 4(6). Both datasets exhibit a proportionality as given by (ojp^z) = — ^d(u/)/dj/ with A « 0.15, 
except for the near-wall region < 25. As a consequence, the average shear Reynolds number 
Re, = \d{uf)/dy\D^/v and the average particle rotation Reynolds number Re^ — {u}p^z)D'^ /v 
are approximately proportional to each other with the proportionality factor A. Incidentally, the 
particle rotation Reynolds number takes values in the interval < Ren < 1.7 (not shown), with 
the maximum occurring at 0.1. 

Velocity covariances for both phases are shown in figure 5. It can be observed that the recorded 
covariance profiles are of very similar values in both simulations, exhibiting only relatively weak 
differences. In absolute terms, the standard deviation values of the streamwise velocity fluctua- 
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Figure 6: Two-point autocorrelations of fluid velocity fluctuations for (a) streamwise separations 
and (5) spanwise separations. The correlation functions are evaluated from 12 instantaneous flow 
fields in the case of Uhlmann (2008) (dashed lines) and from 85 flow fields in the present simulation 
(solid lines). In both cases only data points in the actual fluid domain are taken into account. 

streamwise component [a — 1), wall-normal (a = 2), spanwise (a — 3). The 

data corresponds to a wall-parallel plane at y/h ~ 0.0938 (?/+ « 22). 



tions of both phases ((u^m'^)-^/^ and {u'^u'^^/'^) are reduced on average by 0.007 Mf, and 0.01 u;,, 
respectively. This reduction is consistent with the higher degree of decorrelation of fluid velocity 
data in the extended domain, as discussed below. The curves of the wall-normal and spanwise 
fluid velocity fluctuation intensities in figure 5(a) both nearly collapse with their counterparts from 
the simulation of Uhlmann (2008). On the other hand, the r.m.s. particle velocity fluctuations in 
the wall-normal and spanwise directions are both slightly larger in the present simulation than 
in the previous one (on average by 0.005m;, and 0.01 Ub, respectively). The proflles of the fluid 
Reynolds stress (cf. figure 56) feature a negative peak of slightly smaller amplitude in the present 
dataset, while the values of the covariance between streamwise and wall-normal particle velocity 
fluctuations {u'pv'p) are approximately the same in both datasets in the vicinity of the negative 
peak at y/h ~ 0.15. Further into the core of the channel (i.e. for y/h > 0.4) the present dataset 
exhibits slightly higher negative covariance values, both for (u'j^v'j,) and {u'pv'p). The observed 
differences in the covariance values between the two data-sets is at least in part due to statistical 
uncertainty, as the previous dataset is found to exhibit somewhat noisier profiles (cf. particularly 
the curve for (upVp) in figure 56). However, only when additional simulations with a much longer 
sampling interval are available will it be possible to settle the question about a systematic trend 
in the fluctuation amplitudes. 

Two-point correlations of the fluid velocity fleld evaluated at a wall-distance of y/h = 0.1 
(y+ = 22) are shown in flgure 6. Correlation values for the wall-normal and spanwise velocity 
components are found to be only very little affected by the prolongation of the domain as the 
smaller domain size already yielded a reasonable decorrelation at separations of the order of half the 
fundamental period. On the other hand, fluctuations of the streamwise fluid velocity component 
decorrelate somewhat more rapidly with streamwise separations in the longer box than they do 
in the shorter one. Furthermore, the correlation function Rn in the present simulation reaches 
values close to zero at the largest separations r^/h w 8. Therefore, it makes sense to compute 
streamwise integral length scales {L^aa — ^^''^^ Raa^fx) which take the following values in the 

present case: /h = 0.68, 1^ ^ 0.11, l'"^^ /h = 0.19. Concerning spanwise separations, the 
results from the original and the streamwise-extendcd domain are similar, with the exception of 
Rii which is visibly less noisy in the present data-set. 

Pre-multiplied one-dimensional spectra of the fluid velocity fleld are shown in flgure 7, pro- 
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Figure 7: Premultiplied spectra of the three fluid velocity components as function of the wave- 
length, corresponding to the data in figure 6 (wall-parallel plane at y+ = 22). Each curve is 
normalized to a maximum value of unity in order to emphasize the frequency content, (a) shows 
streamwise spectra, and {h) spanwise spectra. Line-styles and color-coding as in figure 6. Note 
that the particle diameter corresponds to D/h ~ 5 ■ 10^^. 



viding an alternative view of the correlation data of figure 6. In order to compute the spectra, 
velocity values at nodes inside the particles were set to zero, which leads to the well-known 'step- 
noise' at wavelengths around and below the particle scale (Parthasarathy and Facth, 1990), i.e. 
at \/h K, D /h — 0.05. The large- wavelength end of the spectrum, however, is not affected by the 
discontinuities at the phase-interfaces. The spectra are normalized by their respective maximum 
value, in order to enable a comparison of the frequency contents. Concerning the streamwise com- 
ponent of the fluid velocity field, it can be seen in figure 7 that the peak of the energy spectrum 
is captured both in the longer and shorter boxes. However, the decay with wavelength is more 
complete in the larger domain. Nevertheless, the data suggest that still a considerably larger 
streamwise period would be needed in order to allow for a complete energy decay at the large- 
scale end of the spectrum. Unfortunately, a quantitative estimate of the required length cannot 
be deduced from the present data. 

Figures 8 and 9 show streamwise velocity fluctuations from instantaneous flclds of simulations 
in the original domain of Uhlmann (2008) and the extended computational domain of the present 
simulation. Visual inspection of these and other snapshots corroborates the above finding that 
the large-scale streak-like structures (first revealed in Uhlmann, 2008) still exist in the present 
case with a prolonged domain. However, they do not appear to span the entire box length of 
L^/h = 16. 

3.1.2 Lagrangian particle velocity correlations 

Uhlmann (2008) found that Lagrangian velocity correlations were strongly affected by the fiow 
structures with the largest streamwise extension. For the present case, figure 10 shows the auto- 
correlation of particle velocity components along the trajectories. Please note that the correlation 
data is averaged over all particles and integrated over the full observation interval (cf. definition in 
equation 7 of Uhlmann, 2008), which means that no distinction is made concerning the particles' 
instantaneous wall-distance. It had been observed by Uhlmann (2008) that the finite (periodic) 
domain in conjunction with the presence of very long-lived fiow structures leads to a statistical bias 
manifesting itself in the form of successive peaks in the Lagrangian particle velocity autocorrelation 
curves at intervals corresponding to an average return time (i.e. equivalent to the domain length 
divided by the apparent velocity lag u^ei ~ Ub). In the present case we still observe repeated peaks 
in the correlation functions shown in figure 10, but now at intervals of approximately 16 - again 
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Figure 8: Instantaneous three-dimensional 
isosurfaccs of streamwise velocity fluctuations 
u' = Z.&Ur (equivalent to O.Swh). The 
graph (a) corresponds to the simulation in a 
shorter box Uhlmann (2008), (&) is from the 
present simulation in a streamwise-elongated 
box. The view is directed into the wall. Please 
note that only every eighth grid point in each 
direction was used. 



Figure 9: As figure 8, but showing negative- 
valued surfaces at the same magnitude, i.e. 
u' = — 3.6mt-. 
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Figure 10: Lagrangian particle velocity autocorrelation as a function of the separation time t: 

streamwise velocity component (a = 1); - - wall-normal (a = 2); spanwise (a = 3). 

The results previously obtained in a shorter domain (Uhlmann, 2008) are indicated by dashed 
lines. The graph in (6) shows a close-up of the same data as (a) for small separation times. 



consistent with the expected return time {Lx/ub). Furthermore, a faster initial decorrelation of 
the streamwise particle velocity component is observed when using the enlarged box. This finding 
is consistent with the fact that the spatial decorrelation of the corresponding fluid velocity field 
is faster in the longer box (cf. figure 7), although it is still not complete in the present spatial 
domain. 

Despite the statistical bias due to the finite streamwise period, the short-time behavior of the 
auto-correlation in the present case is expected to be well represented. As shown in figure 10(6), 
the initial decay of the correlation function is fastest for the wall-normal component, slightly 
slower for the spanwise component and significantly slower for the streamwise component of the 
particle velocity. The corresponding time scale (Taylor micro-scale, i.e. the intercept with the 
horizontal axis of the osculating parabola at zero separation) measures 0.72 T^, 0.35 Tb and 0.40 Tf, 
for the streamwise, wall-normal and spanwise velocity components, respectively. This result is in 
qualitative agreement with the auto-correlation of fluid particles in single-phase turbulent channel 
flow (Choi, Yeo and Lee, 2004). The explanation for the directional differences put forth by these 
authors is based upon the anisotropy of the near-wall coherent structures (velocity streaks and 
quasi-streamwise vortices) . Despite the presence of particle wakes and additional large scales, the 
flow structure in the present case is similar to single-phase channel flow. Therefore, we expect a 
similar cause to be responsible for the anisotropic autocorrelation of the particle velocity in the 
present case. 

3.2 Voronoi analysis of spatial particle distribution 

A number of techniques have been established for the purpose of characterizing the spatial struc- 
ture of the dispersed phase. In the precursor study (Uhlmann, 2008) conventional box-counting 
(Fessler, Kulick and Eaton, 1994), nearest-neighbor statistics (Kajishima, 2004), as well as genuine 
clustering detection algorithms (Wylic and Koch, 2000; Melheim, 2005) were employed. Based 
upon these measures, it was concluded by Uhlmann (2008) that no significant instantaneous ac- 
cumulation of particles takes place. 

A new technique based upon Voronoi tessellation has recently been proposed in the context of 
particulate flows (Monchaux ct al., 2010; Monchaux, Bourgoin and Cartellier, 2012). Figure 11 
shows an example in two dimensions, where starting with given particle center positions (the 
'sites') the space is covered with cells which have the property that each point inside the cell is 
closer to the cell's site than to any other cell's site. As a consequence, the inverse of a Voronoi 
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Figure 11: (a) Example of a Voronoi tesselation of a bi-periodic (wall-parallel) plane, performed 
with respect to the locations indicated as black dots. Note that the cells are continued periodically 
in order to guarantee a space-filling and non-redundant tesselation. The periodic boundaries are 
indicated by dashed lines, (b) Close-up of a Voronoi tesselation in three-dimensions, as performed 
in the present case. 

cell's volume is an indicator of the local particle concentration. A statistical analysis of Voronoi 
tesselations can be performed by computing the p.d.f. of cell volumes, followed by a comparison 
with reference data for random particle positions (Monchaux et al., 2010, 2012). Compared to 
previous approaches for characterizing the spatial distribution of the dispersed phase, Voronoi 
analysis offers two key advantages: (i) richness of geometrical data, offering the possibility to 
compute various derived diagnostic quantities; (ii) computational efficiency due to the availability 
of a fast algorithm for the tessellation. In the following, we present the results from an application 
of this technique to our present data-set. 

A random distribution of points in unbounded space yields a Gamma distribution for the 
Voronoi volumes (Ferenc and Neda, 2007). However, since our particles are of finite size, a ran- 
dom positioning of particle centers without additional constraints can lead to non-physical overlap 
of particle boundaries. Furthermore, the presence of domain boundaries (walls) can further modify 
the p.d.f. of Voronoi cell volumes. Therefore, we have numerically determined the Voronoi tessel- 
lation of particle positions which have been generated randomly, with the additional constraint 
that no (particle/particle, particle/wall) overlap is obtained. The corresponding p.d.f.'s (which 
are indeed close to a Gamma distribution) are added to the graphs for the purpose of comparison 
with the actual DNS data. 

An additional issue arises in the present context due to the spatial inhomogeneity of the average 
particle concentration. In order to properly account for this variation (i.e. non-constant average 
volume of Voronoi cells) in the wall-normal direction, we have divided the channel height into 
20 intervals of equal width. Statistics of Voronoi cells were then performed individually for each 
'bin', i.e. computed over all Voronoi cells whose site is located inside the corresponding wall-normal 
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Figure 12: Pdf of Voronoi cell volumes. After three-dimensional Voronoi tesselation, pdfs of the 
cell volumes are computed in 20 uniform bins across the channel width. The graph in (a) shows the 
pdf corresponding to the first bin near the wall, (b) is for the bin adjacent to the channel centerline. 
The present DNS data is shown by a solid black curve; the dashed curve corresponds to randomly 
placed (non-overlapping) particles. Data is averaged over 2200 instantaneous snapshots. 



interval. 

A graph of the p.d.f. of Voronoi cell volumes based upon an ensemble of particle positions at 
2200 different times in our present simulation is depicted in figure 12, for the interval adjacent 
to the wall (figure 12a) and for the centerline (figure 126). In studies of particulate flows which 
exhibit significant preferential concentration, an appreciable deviation of the Voronoi cell area 
p.d.f. from the random case is observed (Monchaux ct al., 2010), characterized by larger-than- 
random probabilities of finding very small and very large Voronoi cells while roughly maintaining 
the overall shape of the p.d.f.. Contrarily, our data in figure 12 shows the opposite behavior: very 
small Voronoi cells and very large Voronoi cells are less probable than in the case of randomly 
placed particles. In other words, large deviations from the average cell volume are found to be less 
probable than in a case of randomly placed particles. Therefore, the state of the dispersed phase 
can be characterized as more ordered than random, slightly (but significantly) tending towards 
a homogeneous spatial distribution. It should be noted that the results of Uhlmann (2008) for 
the average distance to the nearest particle (cf. figure 23 therein) did indicate a slight deviation 
from a fully random particle ensemble, tending towards the value for a homogeneous particle array. 
Box-counting and direct cluster identification, however, did not pick up noticeable differences with 
respect to randomness. 

The deviation from a random state as presently observed from Voronoi cell volume statistics 
is consistent with particles having a weak tendency to form a regular pattern. While the extreme 
case of particles forming a uniform cubical lattice would yield a Dirac distribution (with a pulse 
at the average cell volume), the present data shown in figure 12 only weakly deviates from the 
purely random reference curve. On the other hand, preferential particle concentration outside of 
vortical structures can indeed be ruled out in the present case, confirming the previous findings 
of the precursor study (Uhlmann, 2008). 

Additional quantities of interest can be deduced from the Voronoi tessellation. In our case, 
where there exist preferred spatial directions due to gravity and the general non-isotropy of channel 
flow, particle accumulation (if any) is not expected to take place in an isotropic manner. In 
particular, possible elongated particle chains with a preferential orientation can be expected to 
lead to non-isotropic statistics of Voronoi cells. Therefore, we have computed the aspect ratio 
of Voronoi cells dividing the largest streamwise extension by the largest wall-normal ly and 
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Figure 13: Pdf of Voronoi cell aspect ratios, defined as the maximum extension of a given 
cell in the respective coordinate directions: streamwise/wall- normal (A^) are shown in blue, 
streamwise/spanwise {A^ ) are shown in red. As in figure 12, graph (a) shows the pdf corresponding 
to the first bin (out of 20) near the wall, (6) is for the bin adjacent to the channel centerline; solid 
lines correspond to DNS data, dashed lines to random particle positions. 

spanwise extensions Iz of each Voronoi cell, viz. 

Al = j^ Va = 2,3. (2) 

The corresponding p.d.f.'s are shown in figure 13 (again for two different wall-normal intervals). 
Adjacent to the solid surface, the streamwise/spanwise and the streamwise/wall- normal aspect 
ratios are clearly distinct: the latter having a mean of 1.74 due to the constraint by the wall- 
boundary. In the center of the channel the constraint by the wall is no longer felt and the curves 
for the streamwise/spanwise and the streamwise/wall- normal aspect ratios coincide. Except for 
the first interval adjacent to the walls, deviations of the DNS data from the curves for random 
particle arrangements are only observed as an increased probability of very small values Ay and 
AY , i.e. Voronoi cells which are squeezed in the streamwise direction are found more frequently 
than in a random arrangement. This small but systematic difference (cf. figure 135) implies that 
particles have indeed a weak tendency to align along the streamwise direction. It can be speculated 
that such an alignment is induced by the sheltering effect of particle wakes, which is known to 
cause trailing particles to approach leading particles (Wu and Manasseh, 1998; Fortes, Joseph and 
Lundgren, 1987). 

3.3 Particle acceleration statistics 

In this section, we study the Lagrangian acceleration statistics of the solid particles in the present 
flow. In recent years much attention has been given to the Lagrangian acceleration statistics of fluid 
and solid particles, see for example the review by Toschi and Bodcnschatz (2009). Apart from being 
relevant to a number of applications, the growing interest has been mainly due to the appearance of 
new experimental methods and the abundance of direct numerical simulations employing the point- 
particle approach. Most previous studies have been performed for homogeneous flows. Exceptions 
are, for example, the study of Lagrangian acceleration statistics of sub-Kolmogorov size inertial 
particles in a turbulent boundary layer, by means of experiments (Gerashchenko et al., 2008) 
and DNS (Lavezzo, Soldati, Gerashchenko, Warhaft and Collins, 2010). Finite size effects on 
Lagrangian particle acceleration statistics have mainly been studied experimentally, using both 
neutrally buoyant and heavy particles (Qureshi et al., 2007, 2008; Xu and Bodenschatz, 2008; 
Brown et al., 2009). Two contributions based upon DNS studies should be mentioned. Calzavarini 
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Figure 14: Mean particle acceleration, as a 
function of the wall distance. — o— , stream- 
wise component (flp i); — o— , wall-normal com- 
ponent (ap,2)- The lines without symbols 
show the wall-normal gradients of Reynolds 

stresses: , d{u'v') /dy\ - - , d(w'w')/dy. 

Wall units have been used as the normaliza- 
tion scale {aref ~ w^/i^). 




Figure 15: Standard deviation of particle ac- 
celeration as a function of the wall distance. — 

— , streamwise component; , wall-normal 

component; , spanwise component. 



et al. (2009) have extended the point-particle approach by including Faxen corrections, while 
Homann and Bee (20010) have simulated the motion of a single, resolved, finite-size particle, both 
in forced homogeneous-isotropic turbulence. In both cases, however, there was no apparent slip 
velocity due to the exclusion of gravity in the former case and due to matched density in the 
latter. 

Before turning to the present results, it should be pointed out that we have filtered from the 
particle force data all records corresponding to particle collisions (i.e. those where there is a force 
contribution Fc in equation la stemming from the artificial repulsion force). 

Since the present flow is inhomogeneous in the wall-normal direction, we start by considering 
the mean Lagrangian acceleration of the particles, a quantity which is zero in statistically station- 
ary homogeneous flows. Recall that in turbulent channel flow the mean streamwise (wall-normal) 
Lagrangian acceleration of fluid particles turns out to be equal to the wall-normal gradient of 
(u'fv'f) ((^/^/)) (Yeo, Kim and Lee, 2010). Fig. 14 displays the mean streamwise (cpi) and wall- 
normal (ap2) acceleration of the solid particles. Also for comparison the mean streamwise and 
wall-normal Lagrangian fluid particle accelerations (a/i) = dy{u'jv'j) and (a/2) = dy{v'jv'j) are 
shown. The streamwise and wall-normal acceleration components are non-zero near the wall, 
while the mean spanwise acceleration is zero everywhere (not shown). In the central region, (cpi) 
is roughly zero, as can be expected from an equilibrium between buoyancy and drag forces. The 
presence of the wall alters this equilibrium and a negative peak in (opi) appears at y"*" « 10, which 
is significantly larger in magnitude than the corresponding peak observed for (a/i). Furthermore, 
the mean particle acceleration (flpi) presents a milder positive peak at y"*" « 30. Concerning 
the wall-normal acceleration, (ap2) exhibits slightly positive values in the interval 15 < < 40 
(equivalent to a mean lift acting towards the channel center) which are smaller in magnitude than 
those presented by the fluid counterpart (0/2)- Closer to the wall (2/+ < 15) the mean wall- 
normal particle acceleration takes slightly negative values. The peak in (opi) might be explained 
by the following mechanism. Similar to turbulent fluid motion, inertial particle motion near the 
wall exhibits a preference for Q2 and Q4 events (in the terminology of quadrant analysis of the 
streamwise and wall-normal velocity fluctuations, where Q2 refers to ejections of low-speed fluid 
away from the wall and Q4 to an inrush - or sweep - of high-speed fluid), as shown by Uhlmann 
(2008). As a consequence, significant average values of the cross-correlation (i.e. 'Reynolds stress') 
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Figure 16: Mean and rms hydrodynamic particle force components normalized by a reference 
force F^^^y — \pfApU^^l^ (where Ap = nD'^ /A), (a) shows wall-normal profiles, while in (6) the 
data is presented as a function of the local Reynolds number based on an approximation for the 
relative velocity. The red (blue) line is for the rms particle force in the wall-normal (spanwise) 
direction, the relative velocity in the reference force being taken as the rms fluid velocity, i.e. 
"rel ~ ^^/^/)^^^ ("re/ ~ {w'^w'j)^/'^ ). The magcnta-colored line indicates the average particle 
force in the streamwise direction, the relative velocity being taken as the mean apparent slip 
velocity w^.^; = {up) — (uf). The cyan (green) colored lines indicate the mean plus (minus) the 
rms particle force in the streamwise direction, normalized by a reference force based upon the 
apparent slip velocity plus (minus) the rms streamwise fluid velocity. Concerning the graph in 
(6): the horizontal axis represents the Reynolds number based upon the particle diameter and the 
respective relative velocity u';^^^. The black solid line indicates the standard drag law for spheres 
in uniform flow (Clift et al., 1978, table 5.2). The symbols correspond to the data in (a) on the 
channel centerline, using the same color coding. 



(UpVp) arise. Since the particles have larger inertia than a corresponding blob of fluid, it can 
be expected that it takes longer for the particles to adjust to the surrounding fluid conditions, 
and consequently the correlation values can be expected to exceed those of the fluid counterpart, 
as is indeed the case (cf. figure 326). Now, since particles arriving in the near- wall region from 
larger wall-distances carry on average an excess axial velocity value (Q4 events), these particles 
will tend to experience a smaller amount of positive streamwise force (drag) . As a consequence, 
the (negative) submerged weight will exceed the positive drag and a negative particle acceleration 
in the vertical direction will result. The opposite is expected for particles being ejected away 
from the wall (Q2 events). The location where the average acceleration changes sign from negative 
values to positive ones {y/h 0.1, ^ 20) coincides approximately with the location where the 
gradient of the Reynolds shear stress changes sign; it is also an upper bound for the near-wall 
region in which appreciable wall-normal gradients of the wall-normal particle velocity fluctuation 
energy {v'pv'p) exist (cf. figure 32). As a consequence, the turbophoretic effect (Caporaloni et al., 
1975; Recks, 1983; Uhlmann, 2008) can be expected to lead to a preponderance of negative particle 
acceleration for y/h < 0.1, as observed in figure 14. 

Fig. 15 shows the r.m.s. values of the three components of the particle acceleration scaled in wall 
units. The values obtained are comparable in magnitude to the r.m.s. fluid particle acceleration 
in single-phase channel flow (Yeo ct al., 2010) and to the r.m.s. acceleration of inertial particles in 
a turbulent boundary layer (Gerashchenko et al., 2008). The three components present a similar 
trend: a peak near the wall, a dip in the buffer layer and a gradual increase towards a mild 
maximum and finally a plateau in the central region. The r.m.s. of streamwise acceleration is 
larger than the other two components by about 50% except near the wall (y+ < 10), where the 
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Figure 18: Kurtosis K{6u^^) of the Lagrangian velocity increments of solid particle velocities 
5up i (j) = i + ''") ^ i (0 ' plotted as a function of the separation time r . The three coordinate 
directions are color-coded as in figure 15. For each component the dashed line indicates the 
asymptotic value K{u'p^^)/2 + 3/2. 



r.m.s. of the wall-normal component presents maximum values. 

Further insight can be obtained by studying the r.m.s. of the hydrodynamic forces on the 
particles, since by Newton's law studying the force is equivalent (up to a constant coefficient) 
to studying the acceleration. We distinguish here the streamwise direction and the transverse 
directions. For the streamwise direction the mean drag force is approximately balanced by the 
gravity force, while for the transverse directions the mean lift and spanwise forces are zero, except 
for the lift very near the wall, as discussed above. We consider, therefore, for the transverse 
directions the r.m.s. of the hydrodynamic force and for the streamwise direction we consider first, 
the mean hydrodynamic force, and second the mean hydrodynamic force plus (or minus) the r.m.s. 
of the hydrodynamic force. For each direction, we use as a reference force F^''^^ — ^PfApU^^lf 

where Ap = ttD^/I and u"^li is a reference velocity characteristic for the relative fiow in the 
ith coordinate direction. We estimate the relative velocity scale as follows; for the transverse 
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Figure 19: Lagrangian autocorrelation of hydrodynamic force acting upon the particles, plotted 

as a function of the separation time t: , streamwise (a = 1), - - , wall-normal (a = 2), 

, spanwise (a = 3) components. The graph in (6) shows a close-up of the same data as (a) 

for small separation times, including scaling of the abscissa in bulk time units and in wall units. 

directions the r.m.s. fluid velocity is employed, viz. u^Ji = {v'^v'j)^^'^ and m^.^] = {w'^w'j)^/"^; in 
the streamwise direction, the relative velocity is taken as the mean apparent slip velocity (plus 
or minus the r.m.s. fluid velocity), uj^^] — uiag (or — uiag ± (mjU^)^/^), according to whether 
the mean or the mean plus/minus the r.m.s. value of the hydrodynamic force are considered. 
Using this normalization, the mean and r.m.s. hydrodynamic forces are shown in Fig. 16(a). All 
curves are rather flat for y/h > 0.2, and they can be interpreted as force coefficients. Defining 

(i) ■ . . . . . (i) 

a particle Reynolds number RejJ in each direction based on the respective relative velocity u^^j, 

(i) 

and plotting the force coefficients from the central region of the channel as a function of Re)~i , as 
shown in Fig. 16(6), it turns out that the force coefficients are consistent with the standard drag 
law for spheres in uniform flow (Clift ct al., 1978, table 5.2). We have checked that the scaling 
used for the transverse direction leads in the case of the streamwise direction to a force coefficient 
which is not consistent with the standard drag law. As might be expected, the interpretation 
of the normalized hydrodynamic force in terms of a drag coefficient (as shown in Fig. 166) does 
provide an explanation for the large differences in r.m.s. aceleration/force observed when simply 
normalizing all components in wall units (cf. Fig. 15). However, given the rather crude assumptions 
it is indeed remarkable that the mean streamwise force coefficient in Fig. 166 turns out relatively 
close to the value given by the standard drag law valid for a fixed sphere in uniform fiow. 

Next we consider the p.d.f.'s of Lagrangian particle velocity increments i{T) = u'^i{t + 
t) — Mpj(i). Fig. 17 displays the p.d.f.'s of streamwise and wall-normal increments for various 
time lags. Similar graphs have been reported by Mordant ct al. (2002) and Qurcshi ct al. (2007), 
among others. The figure shows that for increasing time increments the p.d.f.'s become more 
like Gaussian distributions. For short time increments, of the order of the viscous time scale, 
pronounced tails are present, which are the signature of intermittent Lagrangian dynamics. It 
is interesting to note the skewness that appears in the p.d.f. corresponding to the streamwise 
velocity increments, a feature that has not been observed in previous investigations in other fiow 
geometries (Mordant et al., 2002; Qurcshi ct al., 2007). This point will be further discussed below. 
The decrease of the tails of the velocity increment p.d.f.'s with increasing separation time r can 
be gauged by examining the temporal evolution of the kurtosis of 5^^ ;(''"), as shown in figure 18. 
It can be seen that the initial decrease is fastest for the streamwise direction, while the curves 
for the two horizontal directions are similar to each other, but with a slight shift in the value 
of the kurtosis. The three components exhibit slight 'bumps' at multiples of \&Tb corresponding 
to the finite-box-size bias (cf. discussion in § 3.1.2), the curves reaching asymptotic values which 
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Figure 20: Normalized probability density function of hydrodynamic forces acting on the parti- 
cles, averaged over the full channel height, , streamwise force component; , wall-normal 

component; , spanwise component; , log-normal fit proposed in Qureshi et al. (2008). 




Figure 21: Higher-order moments of particle acceleration: (a) skewness S, (6) kurtosis K. Lines 
as in figure 15. 



are close to Gaussian for the streamwise direction (K{Sup(T = AOTb) = 3.0) and somewhat larger 
for the two horizontal directions {K{6y^{T = 40Tb) = 4.2 and KiS^opir = •iOTt) = 3.6). It can 
be shown that if the signals at time t and t + t are completely decorrelated, the kurtosis of the 
increment will have a value of K{Up ^)/2 + 3/2 (note that this limit is indeed approached by our 
data at long times as can be seen in figure 18). As a consequence, the different long-time limits 
of K{du i) reflect the differences with respect to intermittency of the three components of the 
particle velocity (cf. the one-time p.d.f.'s shown in flgure 17 of Uhlmann, 2008). 

Let us now turn to the analysis of correlation times of the Lagrangian statistics. Fig. 19 shows 
the Lagrangian autocorrelation of the hydrodynamic forces acting on the particles (which is equiv- 
alent to the autocorrelation of particle acceleration). It can be observed that the autocorrelation 
functions first cross the zero-axis after several viscous time units (the streamwise component after 
Iv/u^, the wall- normal and spanwise components after approx. 3j^/u^), then taking small nega- 
tive values which return to zero on a longer time scale of the order of 2T(,. The initial decay can 
be characterized through the Taylor micro-scale, which is found to measure 2.95, 0.92 and 0.80 



19 



viscous time units for the streamwise, wall-normal and spanwise force components, respectively. 
The figure shows that although the initial decorrelation strongly differs among the three force 
components, their long-time behavior (for separation times r/Tf, > 0.7) is roughly equivalent. A 
very similar directional dependence of the Lagrangian autocorrelation at short separation times 
has also been observed for the acceleration of fluid particles in single-phase channel flow (Choi 
et al., 2004). As in the case of Lagrangian velocity autocorrelations (discussed in § 3.1.2), these 
authors attribute the effect to the anisotropy of the near-wall coherent structures. In view of the 
striking similarity between the present data (figure 196) and the fluid particle data (figure 13, 
Choi ct al., 2004), we believe that the presently observed difference in the initial decorrelation 
time scale can be attributed to the anisotropic nature of the wall-bounded turbulent flow scales. 

In the limit of very short time increments, the p.d.f.'s of Lagrangian particle velocity incre- 
ments, discussed above, become the p.d.f.'s of Lagrangian particle acceleration. Qureshi et al. 
(2008) have shown that acceleration statistics of finite size inertial particles are found very robust 
to size and density variations, the influence of which is mostly carried by the acceleration variance. 
They found that the shape of the normalized p.d.f. remained unchanged over the whole range of 
sizes and density ratios explored. The p.d.f.'s shown in Fig. 20 conflrm this result in the present 
case for the transverse directions (wall- normal and spanwise). Both curves follow the lognormal 
fit proposed by Qureshi et al. (2008) relatively closely over several decades. On the other hand, 
the streamwise acceleration deviates from this curve, showing significantly higher probabilities for 
positive (upward) acceleration fluctuations. It is not clear why this p.d.f. is positively skewed. One 
possible mechanism that provides positive skewness is as follows. Assume that a non-linear drag 
law holds instantaneously, i.e. Cd — f{Reinst), where Re.mst is based on the instantaneous relative 
velocity felt by the particle, and that the drag coefficient decays more slowly with Reynolds than 
^^7nsf such a case, positive velocity fluctuations would lead to a smaller decrease in the drag 
coefficient than the corresponding increase due to negative velocity fluctuations of the same mag- 
nitude. As a consequence the drag force (as well as the corresponding acceleration) can become 
positively skewed even when the velocity fluctuations are symmetric w.r.t. the mean. The data 
in Fig. 16(a) is consistent with this model: it can be seen that the streamwise force coefficient 
obtained using the mean drag plus the r.m.s. drag is closer to the force coefficient of the mean 
drag than the force coefficient obtained using the mean drag minus the r.m.s. drag. We have fur- 
ther tested the hypothesis based upon non-linear drag by using a Gaussian relative velocity p.d.f. 
as input to the standard drag law. Although this simple model indeed yields positively skewed 
hydrodynamic forces (figure omitted), the magnitude of the skewness obtained is smaller than the 
one obtained in the simulation. Therefore, this issue deserves further investigation. 

The force/acceleration p.d.f.'s shown in Fig. 20 have been computed without taking into 
account the inhomogeneous nature of the flow. In order to illustrate how the wall affects the p.d.f.'s 
we have also computed them in 160 wall-normal slabs of uniform thickness. It suffices to discuss 
the wall-normal variation of the p.d.f.'s in terms of their higher order moments, of which skewness 
and kurtosis are shown in Fig. 21. Only the streamwise component presents significant positive 
skewness which varies from S ~ 0.8 close to the wall to S* « 1.3 in the central region. The wall- 
normal and spanwise force/acceleration components are roughly symmetrically distributed, with 
the exception of the region near the wall where a small positive skewness is obtained. Concerning 
the kurtosis, all three components tend to present higher values in the central region of the channel 
where the streamwise component is larger than the other two components by approximately 60%. 
This indicates that the tails of the p.d.f. of the streamwise component are broader as can be 
appreciated in Fig. 20. The kurtosis of the wall-normal and spanwise components are equal in 
the central region of the channel, indicating that in this zone these two directions are equivalent 
as mentioned above. When approaching the wall, they deviate from each other. The kurtosis of 
the spanwise component remains more or less constant at a value of if « 6 up to the immediate 
near- wall region where it decreases mildly. The kurtosis of the wall-normal component, on the 
other hand, presents a remarkable peak in the buffer layer (with K « 12) which is even higher 
than the kurtosis of the streamwise component in this region. Closer to the wall the kurtosis of 
the wall-normal component tends towards a Gaussian value. 
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Figure 22: Relative turbulence intensity as a function of the wall-distance: — o — , Ir 



3.4 Particle-conditioned averaging 

In this section we turn our attention to the statistics of the flow in the neighborhood of the 
particles. For this purpose we have carried out averaging of the flow field conditioned upon the 
presence of particles. Before presenting a brief description of the literature on this topic (§ 3.4.1) 
as well as explaining the averaging procedure (§ 3.4.2), let us introduce the parameters which are 
believed to characterize the local flow field. 

In addition to the Reynolds number based upon the average relative flow velocity and the 
sphere diameter (Reiag), the most prominent parameter describing turbulent flow around spheres 
is the relative turbulence intensity, i.e. the ratio between the intensity of the incoming fluid flow 
fluctuations and the apparent slip velocity. Whereas in homogeneous flows the deflnition of a 
relative turbulence intensity Ir is often based upon the three-component turbulence intensity, viz. 
Ir = {{u[u[) /'iY^'^ /uiagi in cases with unidirectional mean flow (in the x-coordinate direction) the 
following definition is commonly employed: Ir = {u' u')^/"^ /uiag- In order to facilitate a relation 
of our current results to previous and future work, wall-normal profiles of the relative turbulence 
intensity according to both definitions are shown in figure 22. These profiles are found to be 
practically uniform across the channel, with values of Ir ~ 0.14 and Ir ~ 0.23. 

3.4.1 Previous work 

A number of previous studies provide useful reference data for the following discussion of our 
present results. Bagchi and Balachandar (2004) have simulated the flow around a single fixed 
particle in an unbounded domain, swept by homogeneous-isotropic turbulence (at Rex = 164) 
with a superposed uniform mean flow. They considered two values of relative turbulence intensity, 
Ir = 0.1 and 0.25, while varying the particle Reynolds number based upon the average relative 
flow from Reiag — 58 to 610. Although structurally different, homogeneous-isotropic flow is still 
relevant to plane channel flow if the central region of the latter configuration is considered, where a 
far less anisotropic state than in the near-wall region can be found. Therefore, the data of Bagchi 
and Balachandar (2004) will be compared to our present results in the vicinity of the centerplane 
of the channel. 

Amoura et al. (2010) have investigated the flow around a fixed sphere swept by higher in- 
tensity incident turbulence (which was approximately homogeneous-isotropic) with experimental 
techniques. They have considered particle Reynolds numbers based upon the average relative flow 
from Reiag = 100 to 1000, while the relative turbulence intensity was varied from Ir = 0.26 to 
0.45. At the same time, the ratio of particle diameter to integral length scale of the incident 



21 




(c) 



z 
R 



z 
R 



id) 



-5 5 10 15 20 25 30 



R 




y'^^'^h = 0.05, 
y(-'>+ = 11, 
y'^^'^D = 1 

j/^'^)//! 0.15, 

y^''>+ = 33, 
yl'^)/!? = 3 

y("V/i = 0.35, 
y(^^+ = 77, 
^(s)/!? = 7 



?/(*V/i = 0.95, 
= 208, 



Figure 23: Contours of the streamwise component of the average relative velocity u in wall- 
parallel planes through the center of the particles, plotted for slabs at different wall-distances y^^\ 
Contourlines are at uniformly distributed values (0 : .1 : 0.9) times the maximum value in each 
plane. The zero-valued contour is plotted in red color. 



turbulence was significantly larger than in most previous studies. 

Wu and Facth (1994) have performed experimental measurements of the flow around fixed 
spheres on the centerline of turbulent pipe flow {Rei, — 17000 up to 50000, based upon pipe 
radius). The particle Reynolds number based upon the average relative velocity was varied from 
Reiag = 135 to 1560. 

Legendre et al. (2006) have simulated the flow around a stationary particle on the centerline 
of turbulent pipe flow [Ret, — 3000 and Re^ = 200, based upon the pipe radius) by means of 
LES. Their particle has a diameter of £)+ ~ 10.23 and a particle Reynolds number based upon 
the average relative velocity of Reiag — 200, providing a parameter point which is not far removed 
from our present study. The prime difference is their relative turbulence intensity Ir = 0.037, 
which is significantly lower than in the present case. 

Zcng et al. (2010) have simulated turbulent channel flow (Rer ~ 180) with a single fixed 
particle, either located in the buffer layer (?/+ « 18) or at the center of the channel. They also 
varied the size of the particle in the range D'^ — 3.5 to 25, such that a total of six different cases 
were considered. Out of that data-set, their case 2 is directly comparable to our data, featuring a 
particle with diameter = 10.7 placed in the buffer layer. 

Please note that all of the studies mentioned above consider the flow around single particles 
which were fixed in space. Consequently, effects of particle mobility and collective effects were 
absent in those cases. We are not aware of previous work involving particle-conditioned averaging 
of turbulent flow-fields at the scale of the particles. One exception is the experimental work of 
Poelma, Westerweel and Ooms (2007) who did measure the conditionally-sampled flow field around 
mobile particles in decaying grid-generated turbulence. However, the near-field results presented 
in that reference are not sufficiently detailed in order to serve for the purpose of comparison in 
the present context. 
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Figure 24: Average streamwise recirculation length in the wake of the particles, shown as a 
function of the wall-distance. Results for fixed single particles from the literature are indicated 
by the following symbols; particle in buffer layer of turbulent channel flow (case 2 of Zeng 
et al., 2010); particle in buffer layer of laminar channel flow (Zeng ct al., 2010); ■, particle 
in isotropic turbulence with mean relative-flow Reynolds number of Reiag = 114 and relative 
turbulence intensity Ir = 0.25 (case 4 of Bagchi and Balachandar, 2004); ■, particle in uniform 
flow at mean relative-flow Reynolds number of Reiag = 114 (Bagchi and Balachandar, 2004). 

3.4.2 Preliminaries 

In the following we present an analysis of the local relative velocity field in the vicinity of the 
particles. For this purpose we have averaged the Eulerian velocity field in a frame of reference 
attached to the particle centers, summing over particles located in 20 wall-normal slabs (width 
of the slabs equal to 2D) as well as over a number of 97 instantaneous fields (yielding more than 
70000 samples per slab). Before turning to the results, let us first make the averaging process 
more precise. 

The instantaneous velocity difference between the two phases, computed with respect to a 
given particle, is defined as follows: 

uW(x,t)=u/(x,i)-uW(t), (3) 

where the instantaneous coordinate relative to the ith particle's center position x.p \t) is given by 

x«(t)=x-x«(i). (4) 

Furthermore, in (3) u/(x, i) refers to the fluid velocity vector fleld and Up •*(<:) is the ith particle's 
velocity vector at a given time t. The average relative velocity u{i,y^-^^) in the slab centered at a 
wall-distance y^^^ is computed from 

£i(i,2;W) = (uW(i,t))p,t, (5) 

where the operator (■)p,t denotes averaging over all particles (whose center position is located in 
the respective slab) and over time (i.e. over a number of snapshots). Please note that averaging 
is performed over the actual volume instantaneously occupied by the fluid. This means that 
the volume occupied by other particles located in the neighborhood of a particular sphere is 
disregarded in the averaging procedure shown in (5). The reader is referred to Appendix A for a 
precise definition of the above averaging operator. 

We can now decompose the instantaneous value of the relative velocity w.r.t. the ith particle 
(3) into an average contribution (5) and a fluctuation u^*^", viz. 

u«(i,0 = ii(i,2/(^') + u«"(i,i), (6) 
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Figure 25: Average relative streamwise velocity u along the streamwise axis through the center 
of the particles. The velocity u is normalized by the apparent velocity lag uiag = \{uf) — {up)\. 
The solid curves correspond to particles being located at different wall-distances, averaged over 
bins centered at: = 11, = 33, = 77, y^''^+ = 208. 



where the average of the fluctuation vanishes identically (i.e. (u^ ' )p t = 0). From the fluctuations 
Ur^"(x, t) implicitly defined in (6) we can compute second moments; for a component in the a- 
direction the definition reads: 

where no summation is implied for greek indices. Furthermore, we can define velocity fiuctuations 
in the particle-centered frame of reference for each of the phases, viz. 

\if"{ii,t) = u/(x,t) - (u/(i,i))p,t , (8a) 
u«"(i) = uW(i)-(u,(t))p,,, (8b) 

such that xif'" = — Up 

3.4.3 Mean relative velocity field 

Figure 23 shows the streamwise component of the average relative velocity, u, in wall-parallel 
planes through the particle center for averaging slabs at different wall-distances. It can be seen 
that the overall wake pattern in those planes is similar across the channel. However, the extent of 
the recirculation region as well as the deceleration (acceleration) upstream (downstream) of the 
particle are found to depend upon the wall distance. In particular, the region where the average 
streamwise velocity exhibits negative values is somewhat smaller at y^"-*"*" = 11 than in the outer 
region. Likewise, the stagnation region upstream of the particle is more compact at this near-wall 
location (j/^'*-'^ = 11) than at larger wall-distances. It can also be observed in figure 23 that there 
are no significant differences between the wake pattern of particles in the logarithmic layer and 
those located in the outer region or at the center of the channel. 

The length of the recirculation region Lg, defined as the streamwise distance from the rear 
stagnation point to the downstream location where the average relative velocity changes sign, is 
shown in figure 24 as a function of wall-distance. Please note that is measured in the wall- 
parallel plane passing through the particle center; therefore, the value of does not necessarily 
reflect the maximum length of the three-dimensional recirculation region, which might be located 
off-center. The figure shows an approximately constant value of Lf./R « 1.72 for wall-distances 
y/h > 0.2 (2/+ > 45), where R — D/2 is the particle radius. For smaller wall-distances, the 
length of the recirculation zone decreases, down to a value of L^/R = 1.48 in the first slab at 
y^''')/h = 0.05 (?/('')+ = 11). 



24 



UdO 




10° 10' 10^ 



i/R 

Figure 26: The streamwise velocity deficit Udo 
along the streamwise axis through the center 
of the particles (on their downstream side). 
Color coding as in figure 25. The dashed 
straight lines indicate decay rates proportional 
to and x^'^ . 
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Figure 27: As in figure 26, but the veloc- 
ity deficit is averaged over all particles lo- 
cated at y/h > 0.2 {y/D > 4), and the result 
is compensated in order to highlight power- 
law regions, i.e. the quantity Udo x" is plotted, 
where: , n = 1; - - , n = 2. 



Figure 24 also includes data points for single fixed spheres in laminar channel and uniform flow 
(i.e. the laminar results reported by Bagchi and Balachandar, 2004; Zeng et al., 2010) at compa- 
rable average-flow particle Reynolds numbers. Since homogeneous-isotropic inflow was considered 
by Bagchi and Balachandar (2004), we relate their data to our results at the channel center. It 
can be seen that the presently obtained recirculation lengths fall below the reference values in 
laminar flow. This result is consistent with previous observations made with fixed spheres swept 
by turbulent flow: both Bagchi and Balachandar (2004) as well as Zeng et al. (2010) observed 
a reduction of the recirculation length due to background turbulence. Compared to the values 
of those authors (cf. black symbols in flgure 24) our present results show somewhat larger recir- 
culation lengths in turbulent background flow at comparable turbulence intensity. In particular, 
the present recirculation length in the buffer layer is approximately 26% higher than the value 
reported for case 2 of Zeng et al. (2010); the present results in the core of the channel are 6% 
higher than case 4 of Bagchi and Balachandar (2004) . Although the differences are not extremely 
large, they might reflect the difference in the physics between our case and the single flxed particle 
configurations. Apart from possible effects of particle mobility, it should be remembered that 
the present flow field is significantly altered by the presence of particles, i.e. particle wakes are 
prominent features. Consequently, particles experience a turbulent flow field which is structurally 
significantly different from both the canonical channel fiow and homogeneous-isotropic turbulence 
swept over fixed particles (Bagchi and Balachandar, 2004; Zeng et al., 2010). 

Figure 25 shows the average relative velocity on the streamwise axis through the particle center. 
It should be noted that the curves do not quite reach the particle surface due to the interpolation 
during particle-centered averaging, as explained in Appendix A.l. In figure 25 the relative velocity 
u is normalized by the apparent lag w/ag- It can be seen that with increasing distance from the 
particle, the average relative velocity tends to unity for all wall-distances. It is also visible from 
figure 25 that for particles located in the slab centered at y'-*''*' = 11 the approach to unity appears 
slightly faster (both on the upstream and the downstream side) as compared to particles at larger 
wall distances. 

In order to further investigate the recovery of the average velocity in the wake, we define the 
normalized average velocity deficit in the wall-parallel plane passing through the particle center, 
viz. 
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Figure 28: Average velocity deficit in the particle wakes as a function of the normalized spanwise 
coordinate, given at two downstream locations: o, x/R = 10; A, x/R = 20. The color code 
indicates the wall-distance of the averaging slab (as given in figure 25). The dashed line corresponds 
to a Gaussian function as defined in (12). 

where the average incoming relative velocity Moo(y^*-') is defined as the maximum (over x) of 
M((i, 0, 0)-^, y(^)). The velocity deficit along the streamwise axis through the particle center, Udo, 
is then simply defined as follows 

Udo(i,y(^^) = Urf(i,2/(^),0). (10) 

Figure 26 shows the downstream evolution of the normalized velocity deficit Udo in double loga- 
rithmic scale for distances of up to 140i?. As observed in previous studies on turbulent flow around 
fixed single spheres (Wu and Faeth, 1994; Legendre et al., 2006; Amoura et al., 2010), roughly two 
regions can be distinguished: a near-wake zone with a decay rate approximately proportional to 
x~^ , and a far- wake with a decay of approximately x"^ . The velocity deficit for particles in the av- 
eraging slab adjacent to the wall (y^*^"*" = 11) is somewhat smaller for distances up to x/R « 10, 
as already observed above. Otherwise, the curves corresponding to different wall-distances are 
equivalent to within statistical uncertainty. It has been observed by Legendre et al. (2006) that 
the change in slope (from to x~'^) takes place at a downstream location where the velocity 
deficit and the turbulence intensity are of the same order (udQ ~ (w'w')^/^). The same is true 
in the present case: the change in slope occurs at approximately x/R — 25, where the velocity 
deficit measures Udo ~ 0.1, which is indeed a value comparable to the turbulence intensity at the 
corresponding wall- normal locations (cf. figure 5a). 

The extent of the regions where power-law behavior is observed can be deduced from figure 27 
which shows the compensated velocity deficit u^o i" with exponents n = 1 and n = 2. In order to 
further increase the number of available samples, the data of all particles located at y/h > 0.2 (i.e. 
slabs centered at y^*-' /h > 0.25) has been averaged. The figure shows that a decay of the velocity 
deficit in the particle wakes according to x~'^ takes place in a region approximately delimited 
by 40 < x/R < 80. By way of comparison, in relatively low-intensity turbulence (/^ — 0.037) 
Legendre et al. (2006) found the x~^ law to hold for distances above X2 = 50i?, while Amoura 
et al. (2010) obtain the same evolution for X2 > 5i? at higher turbulence intensity (/^ > 0.26). 
The present results can therefore be qualified as consistent with the trend exhibited in these two 
studies. 

Let us now turn to the question of self-similarity of the particle wakes. Wu and Faeth (1994) 
and subsequent authors (Bagchi and Balachandar, 2004; Legendre et al., 2006) observed that 
cross-stream profiles of the velocity deficit Ud in the wake of isolated fixed spheres in turbulent 
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Figure 29: Spanwise half- width of the average particle wakes as a function of downstream distance, 
(a) shows data for individual wall-normal slabs, with color-coding as in figure 25. (6) shows data 
averaged over the core of the channel (for ylh> 0.2). The dashed line in (&) is given by the 
function h^jR = 0.04 i/i? -|- 0.65. 



surroundings follow a Gaussian function, viz. 

^-exp(--l^). (11) 

An appropriate length scale in (11) is the half- width h^z defined as the lateral position where the 
following relation holds: 

Ud{z = h^z)/udo = exp(-l/2) . (12) 

Figure 28 shows that the cross-wake profiles of the velocity deficit in the present case do indeed 
follow the Gaussian function (11) reasonably well; this is true for all wall-distances and over a 
substantial axial distance downstream of the particles. The streamwise evolution of the spanwise 
half-width h^z as defined in (12) is shown in figure 29. The graph in figure 29(a) demonstrates 
that the particles' wall-distance has a minor effect upon the wake half-width, as all curves have 
very similar evolutions; again only the averaging slab adjacent to the wall (y^*^"^ = 11) represents 
a small exception, with a computed half-width which is systematically at slightly smaller values. 

In figure 29(6) the core-averaged half-width (i.e. found upon averaging over all particles located 
at y/h > 0.2) is shown. A region of linear growth, i.e. 

h^z/R^ ax/R + P , (13) 

which extends over 15 < x/R < 50 is observed. Let us recall that a linear expansion hzw ~ i is 
different from the case of spheres in uniform fiow, in which laminar wakes obey hzw 5;^^^ and 
turbulent wakes exhibit hzw ~ x^^^- The streamwise evolution in the present case (many mobile 
particles in turbulent background fiow) is quite accurately represented by the coefficient values 
a = 0.04 and f3 = 0.65 in the linear expansion law (13). At much lower relative turbulence intensity 
{Ir = 0.037), Lcgendrc ct al. (2006) have found a value of a = 0.024 for the wake expansion rate. 
Bagchi and Balachandar (2004), on the other hand, have obtained a value of a = 0.135 for their 
cases with relative turbulence intensity of — 0.1, independently of the particle Reynolds number. 
The analytic results of Barnes et al. (2011) for wake spreading in turbulent surroundings do suggest 
a value of the expansion factor (in the linear spreading regime) comparable to the value of the 
relative turbulence intensity, viz. a « 7^. This is apparently not the case in our present flow. 
However, it should be kept in mind that in the core of the channel, the dominant contribution 
of turbulent kinetic energy stems from the particle wakes themselves and not from 'incoming' 
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Figure 30: Energy of covariances between relative velocity fluctuations u'^iu'^i downstream of the 
particles, normalized by the reference kinetic energy given in (16) The color code indicates the 
wall-distance of the averaging slab (as given in figure 25). The dashed line corresponds to the 
data averaged over the core of the channel (for y/h > 0.2). 

turbulence in the classical sense. Therefore, some of the assumptions made in the derivation of 
their model (in particular the homogeneous-isotropic structure of the turbulent flow field) do not 
appear to apply. 

3.4.4 Fluctuations 

By using the definition of the fluctuating velocities of each phase (8) the definition of the covari- 
ances between fluctuating relative velocities (7) can be rewritten as follows: 

u^Jic,y^^^) = (4Jx,t)4Jx,i))p,t + W)p,t - 2(4„(i,i)«W^'(0)p,t . (14) 

This relation shows that the covariances of fluctuations w.r.t. the average defined in (5) result from 
three contributions: (i) covariances of the fluctuating fluid velocity field (conditioned on particle 
presence), (ii) covariances of particle velocity fluctuations, and (iii) covariances between particle 
velocity fluctuations and particle-conditioned fluid velocity fluctuations. For fixed particles the 
last two contributions in (14) vanish identically. Furthermore, for large distances from the particle 
(i.e. for large values of x and/or z) we expect the last contribution (covariances between fluid 
and particle velocity fluctuations) to vanish. In that long-distance limit, the contribution from the 
particle-conditioned fluid velocity fluctuation covariances (first term on the r.h.s. of 14) is expected 
to tend towards the simple (unconditioned) plane-and-time-averaged fluid covariance {u'j ^u'j ^). 
Therefore, we define the following velocity scale: 

which is expected to allow for a reasonable normalization of the covariances of particle-conditioned 
relative velocities (7). From (15) we can derive a reference fiuctuation energy k^'^-^ , viz. 

k7'=<:/<://2. (16) 

The downstream evolution of fiuctuation energy u"jU"j/2 on the axis through the particle center 
is shown in figure 30. At all wall-distances, the fluctuation energy increases from zero at the 
rear particle surface and tends towards unity at large downstream distances. For intermediate 
distances, two regimes can be distinguished. In the near- wake {x/R < 10, corresponding to 
the region where Udo ~ x~^) a rapid increase of fluctuation energy is observed, which evolves 
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approximately linearly with x. Further downstream, the approach towards unity is much slower, 
approximately following a logarithmic dependency on the downstream distance x. 

4 Conclusion 

In the present study we have revisited the case of vertical plane channel flow seeded with finite-size 
heavy particles already investigated by Uhlmann (2008). An additional DNS has been carried out 
with identical parameters as in the main simulation of that reference, except for the streamwise 
period of the computational domain which has been doubled. 

The new dataset has first been compared to the previous simulation data of Uhlmann (2008), 
in order to determine the infiuence of the box size on the largest flow structures. It was observed 
that the columnar flow structures, which are induced by the presence of particles, are still not 
fully decorrelated in the prolonged domain. As a consequence, the Lagrangian auto-correlation 
of particle velocity is still biased by the fact that particles re-encounter long-lived flow structures 
after one return time (based upon the apparent slip velocity and the domain size). Results for 
average Eulerian quantities such as the mean fluid and particle velocity profiles are not strongly 
affected by the streamwise extension of the domain. 

Furthermore, we have analyzed the new dataset with respect to additional aspects not previ- 
ously considered in the context of vertical particulate channel fiow. First, we have conducted an 
analysis of the spatial distribution of the disperse phase based upon Voronoi tessellation. It was 
found that the particles are less disorderly distributed than in a random case, slightly tending 
towards a homogeneous distribution. In addition, examining the aspect ratio of Voronoi cells has 
revealed that the particles exhibit a weak tendency to align in the streamwise direction, a trend 
which may be related to wake sheltering. 

Second, we have carried out an analysis of the statistics related to particle acceleration. Near 
the wall, mean particle acceleration is found to significantly deviate from the mean acceleration 
of fluid particles. The standard deviation of particle acceleration, when normalized in wall units, 
differs among the three spatial components. When expressed as fluctuations of force coefficients 
(with approximate velocity scales based on the mean and rms fluid velocity of the corresponding 
components), this variability can be interpreted in terms of the standard drag law. Concerning 
the temporal correlation of particle acceleration data, we have observed a first zero-crossing after 
several viscous time units, followed by a slower decorrelation over several bulk time units. The 
p.d.f.'s of particle acceleration in the wall-normal and spanwise directions are found to be consistent 
with a lognormal distribution as proposed by Qureshi ct al. (2008). The streamwise component, 
however, deviates from that lognormal fit, exhibiting significant positive skewness. One possible 
explanation for this positively-skewed acceleration p.d.f. is through a non-linear drag mechanism 
which would solely affect the streamwise component since it is the only spatial component with a 
finite apparent slip velocity. This point certainly merits further investigation. 

Finally, we have performed particle-conditioned averaging of the flow field in the vicinity of 
the particles. It was found that the characteristics of the particle wakes in the present case are 
nearly independent of the wall distance, except for the near-wall region {y/h < 0.2). The average 
length of the recirculation zone in the wake of the present particles (with relative flow Reynolds 
number based upon apparent slip velocity of 132 in the core of the channel) is consistent with 
data for single fixed spheres investigated by Bagchi and Balachandar (2004) and Zeng et al. (2010) 
at comparable Reynolds number and turbulence intensity. In particular, our results in the buffer 
layer (compared to case 2 of Zeng et al., 2010) and on the centerline (compared to case 4 of 
Bagchi and Balachandar, 2004) feature a somewhat larger mean wake length. The streamwise 
velocity deficit on the axis through the particles is found to decay as a;~^ in the near- wake and 
as for distances beyond approximately 40 particle radii. As in previous studies involving 
fixed spheres in turbulent fiow (Wu and Faeth, 1994; Bagchi and Balachandar, 2004; Legendre 
et al., 2006), the present average streamwise velocity profiles in the particles' wakes can be fit 
with reasonable agreement to an exponential function. The wake half-width resulting from the fit 
evolves almost linearly with the downstream distance over a considerable interval. The energy of 
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the velocity fluctuations with respect to their particle-conditioned average, which is by definition 
zero on the particle surface, is found to evolve approximately linearly with downstream distance 
in the near-wake, and according to a logarithmic law in the far-wake. 
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Appendix A Averaging procedures 



A.l Particle-centered averaging 

Let us define an indicator function (t)l^.-li{y) which signals whether a given wall-normal position y 
is located inside or outside a particular wall-normal slab with index j, viz. 

"PunW) I else ' ^ ' 

where N^in is the number of slabs used to span the channel width. Similarly, we define an indicator 
function (f)f{x.,t) for the fluid phase: 

^/(^'^)^| else ' (1^) 

where ^f(t) is the part of the computational domain f2 which is occupied by fluid at time t. We 
can now define a discrete counter field n|^], which holds the number of samples obtained through 
averaging at a given grid node with indices i, j, k for a given y-slab with index s, viz. 

m=l 1=1 

In equation (19) the symbol t™ indicates the time corresponding to the mth snapshot in the 
database (comprising a total of N^nap snapshots) and ^fjkit"^) is a coordinate relative to the Zth 
particle's center position at time t™, as defined in (4). 

The actual averaging can now be performed analogously to (19), including a division by the 
local number of samples. For a vector field ^(x, t) we define the averaging operator as follows: 

N,„ap JVp 

^ijk m=l 1=1 

In practice the coordinates x^^^j,(i™) of the grid nodes covering the predefined averaging volume do 
not coincide with the grid used in the direct numerical simulation. This means that the vector field 
to be averaged (^ in equation 20) is not directly available at the coordinates x|^--*j,(t™). Therefore, 
the averaging process defined in (20) involves spatial interpolation, which has been realized with a 
tri-linear formula. As a consequence of the use of the fluid indicator function, the particle-centered 
average fields do not quite reach the particle surface. 

The number of (uniformly spaced) bins for evaluating averages defined by (20) was chosen 
as Nbin — 20. The number of snapshots amounts to Nsnap — 97, yielding approximately 80000 
samples per slab. Further averaging over all bins with y^^^ /h > 0.25 (also referred to as "core- 
averaging" in the main text) yields approximately 650000 samples for the quantities shown in 
figures 27 and 29(6). 



A. 2 Wall-parallel plane and time averaging 

The standard averaging over wall-parallel planes and time (not conditioned upon particle presence) 
can now conveniently be defined using the notation and the indicator functions already introduced 
in A.l. Let us first define a counter of fluid sample points in a wall-parallel plane at a given wall- 
distance with index j, viz. 

E EE'^/(^^^'='*'")- (21) 

m— 1 i=l k=l 

Note that the number of data-sets Nume used for this quantity is much larger than Nsnap in the 
case of particle-centered averaging (equations 19 and 20), since standard averaging is performed 
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Figure 31: (a) Mean velocity profiles ol both phases in the present simulation: , fluid phase; 

• , solid phase. (6) Apparent slip velocity between the phases. The fluid phase data is averaged 
according to the operator defined in (22). 



at runtime, while the latter is carried out at a post-processing stage. Averaging over wall-parallel 
planes and in time, while considering only grid points being located in the fluid domain, is deflned 
as follows: 

^ m—l i — 1 k—1 

Consequently, is a function of wall-distance alone, whereas {C)p,t is a three-dimensional field 
for each wall- normal slab y^^\ 

Note that purely for the purpose of strict comparison between the present simulation results 
and data from Uhlmann (2008), the following averaging operator which does not distinguishes 
between the solid and fluid phases is used ('composite' averaging, cf. figures 3 and 5): 

= nt-ww ^ EE^(-^.-'*'")- (23) 

JVt^meJV^iV^ m=l t=l k = l 

A. 3 Binned averages over particle-related quantities 

Concerning Lagrangian quantities, we employ averages over wall-normal bins, using the indicator 
function given in (17). The sample counter for each bin with index s is computed over the number 
^time available particle fields and summing over all particles, viz. 

f^''^- E E^^»(4'H^'")). (24) 

m=l 1=1 

The binned average (over time and the number of particles) of a Lagrangian quantity is defined 
as follows: 

{Cp){i'^) = ^ T.^Y:<i>'-liylHn)c/\n- (25) 

m=l 1=1 

In this manuscript we have chosen the number of (uniformly spaced) bins for evaluating averages 
defined by (25) as Nhin = 160. The present data corresponds to a total number of approximately 
1.8 • lO'' samples. 
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Figure 32: (a) R.m.s. of velocity fluctuations of both phases in the present simulation, with lines 
corresponding to the fluid phase, and symbols to the particulate phase. The color coding indicates 
streamwise (black), wall-normal (red), spanwise (blue) components. (6) Reynolds shear stress of 
fluid velocity fluctuations as well as corresponding velocity correlation of the particle motion. All 
quantities are normalized in bulk units. The fluid phase data is averaged according to the operator 
defined in (22). 

Appendix B Eulerian statistics for the present simulation 

Figures 31 and 32 show the same quantities as presented in figures 3 and 5 (for the present 
simulation), but defined with the correct averaging operator which only considers grid nodes 
being instantaneously located in the fluid domain (cf. definition in equation 22). These quantities 
are presented for future reference. 
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